home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / GAUSSJ.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  2KB  |  83 lines
  1. PROCEDURE gaussj(VAR a: glnpbynp; n,np: integer;
  2.        VAR b: glnpbymp; m,mp: integer);
  3. (* Programs using GAUSSJ must define the types
  4. TYPE
  5.    glnpbynp = ARRAY [1..np,1..np] OF real;
  6.    glnpbymp = ARRAY [1..np,1..mp] OF real;
  7.    glnp = ARRAY [1..np] OF integer;
  8. in the main routine. *)
  9. VAR
  10.    big,dum,pivinv: real;
  11.    i,icol,irow,j,k,l,ll: integer;
  12.    indxc,indxr,ipiv: glnp;
  13. BEGIN
  14.    FOR j := 1 TO n DO BEGIN
  15.       ipiv[j] := 0
  16.    END;
  17.    FOR i := 1 TO n DO BEGIN
  18.       big := 0.0;
  19.       FOR j := 1 TO n DO BEGIN
  20.          IF (ipiv[j] <> 1) THEN BEGIN
  21.             FOR k := 1 TO n DO BEGIN
  22.                IF (ipiv[k] = 0) THEN BEGIN
  23.                   IF (abs(a[j,k]) >= big) THEN BEGIN
  24.                      big := abs(a[j,k]);
  25.                      irow := j;
  26.                      icol := k
  27.                   END
  28.                END ELSE IF (ipiv[k] > 1) THEN BEGIN
  29.                   writeln('pause 1 in GAUSSJ - singular matrix'); readln
  30.                END
  31.             END
  32.          END
  33.       END;
  34.       ipiv[icol] := ipiv[icol]+1;
  35.       IF (irow <> icol) THEN BEGIN
  36.          FOR l := 1 TO n DO BEGIN
  37.             dum := a[irow,l];
  38.             a[irow,l] := a[icol,l];
  39.             a[icol,l] := dum
  40.          END;
  41.          FOR l := 1 TO m DO BEGIN
  42.             dum := b[irow,l];
  43.             b[irow,l] := b[icol,l];
  44.             b[icol,l] := dum
  45.          END
  46.       END;
  47.       indxr[i] := irow;
  48.       indxc[i] := icol;
  49.       IF (a[icol,icol] = 0.0) THEN BEGIN
  50.          writeln('pause 2 in GAUSSJ - singular matrix'); readln
  51.       END;
  52.       pivinv := 1.0/a[icol,icol];
  53.       a[icol,icol] := 1.0;
  54.       FOR l := 1 TO n DO BEGIN
  55.          a[icol,l] := a[icol,l]*pivinv
  56.       END;
  57.       FOR l := 1 TO m DO BEGIN
  58.          b[icol,l] := b[icol,l]*pivinv
  59.       END;
  60.       FOR ll := 1 TO n DO BEGIN
  61.          IF (ll <> icol) THEN BEGIN
  62.             dum := a[ll,icol];
  63.             a[ll,icol] := 0.0;
  64.             FOR l := 1 TO n DO BEGIN
  65.                a[ll,l] := a[ll,l]-a[icol,l]*dum
  66.             END;
  67.             FOR l := 1 TO m DO BEGIN
  68.                b[ll,l] := b[ll,l]-b[icol,l]*dum
  69.             END
  70.          END
  71.       END
  72.    END;
  73.    FOR l := n DOWNTO 1 DO BEGIN
  74.       IF (indxr[l] <> indxc[l]) THEN BEGIN
  75.          FOR k := 1 TO n DO BEGIN
  76.             dum := a[k,indxr[l]];
  77.             a[k,indxr[l]] := a[k,indxc[l]];
  78.             a[k,indxc[l]] := dum
  79.          END
  80.       END
  81.    END
  82. END;
  83.